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qq ' We investigate the large- iV phase transition of lattice SU(JV) gauge theories in the Wilson formulation, by per- 

, forming a Monte Carlo simulation of the twisted Eguchi-Kawai model. A variant of the multicanonical algorithm 
^\ ■ allows a detailed exploration of the phase transition and a precise determination of the transition temperature. 
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1. Introduction 

The large- AT limit of SU(7V) gauge theories is of 
considerable interest from the phenomcnological 
point of view, and is one of our sources of under- 
standing of non-perturbative QCD. The theoret- 
ical side is no less interesting, allowing e.g. to es- 
tablish precise relationships between Yang-Mills 
theories and string theories. 

A crucial property of large- N field theories is 
factorization: connected Green's functions of in- 
variant quantities are suppressed with respect to 
the corresponding disconnected parts by powers 
of 1/N. One of the most notable consequences 
of factorization is the possibility of constructing 
reduced models, i.e. single-site models which re- 
produce a lattice gauge theory in the N — > oo 
limit @,|. 

Despite the considerable simplifications occur- 
ring for large N, many interesting models have 
not been solved analytically; therefore we must 
resort to approximate techniques, such as numer- 
ical simulations. 

2. The TEK model 

The most promising reduced version of the 
SU(A r ) lattice gauge theory (in the Wilson formu- 
lation) is the twisted Eguchi-Kawai model (TEK) 
H). Let us define a set of 4 traceless SU(iV) ma- 
trices T M obeying 't Hooft algebra 
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where n^ v is an antisymmetric integer- valued ten- 
sor, are the matrices implementing the trans- 
lations by one lattice spacing in the /i direction. 



The observables of the reduced model are ob- 
tained by applying the reduction prescription 



U^x) -» T(x)U„T( X y 



to the corresponding quantity of the full model. 
The action of the TEK model is obtained by re- 
duction of the Wilson action: 

I3N 2 E(U) = S T ek(U) = 

N/3 Tr [Z^UnVvUpl + h.c] . 

For a proper choice of T M , the Schwinger-Dyson 
equations of the TEK model reproduce in the 
large-A^ limit the Schwinger-Dyson equations of 
the Wilson lattice gauge theory. We adopt the 
choice of Ref . : N is constrained to be a per- 
fect square, N = L 2 , and n Ufi = L for all v > fi. L 
takes the role of the lattice size, and the number 
of degrees of freedom is proportional to L 4 . 

The TEK model develops a first-order phase 
transition for N — > oo. It corresponds to the 
large- N limit of the first-order phase transition 
of SU(iV) lattice gauge theories. We will study 
the phase transition of the TEK model by Monte 
Carlo simulation. 

3. The multicanonical algorithm 

The choice of updating algorithm is extremely 
important in the neighborhood of a phase tran- 
sition. Most numerical work on the TEK model 
adopted the heat-bath updating algorithm (HB) 
devised in Ref. j^j. The HB algorithm is gener- 
ally quite efficient, but near the phase transition 
it is plagued by the familiar supercritical slowing 
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down; according to the experience accumulated 
in the last decade, this problem can be solved by 
devising a suitable multicanonical algorithm jfj). 

Since we are studying a temperature-driven 
phase transition, it is natural to choose a re- 
weighting function depending only on the energy: 
we generate configurations according to the prob- 
ability distribution 

P{j3,U) oc w(E(U)) exp(-f3N 2 E(U)) 

and compute the canonical expectation value of 
an observable as 

(0)=J2^ 1 (E(c))0(c) /5>- l (25(c)). 

c c 

If we can find a is reasonably accurate ansatz p 
to the energy distribution p(E), we can construct 
the reweighting function 

w{E)=\/p(EJ), E<E_, 

w(E) = l/p(E), E_<E<E + , 

w(E) = l/p(E+), E+ < E, (1) 

which flattens the probability distribution be- 
tween the two peaks, and we can obtain a reason- 
able tunneling rate between the two vacua. The 
resulting multicanonical distribution can be sim- 
ulated using an efficient Metropolis procedure. 
Our first ansatz was 

Pbl(E) = a+g+(E) + a-g-(E) +j8{E-E+) 

x9(E_-E) (1 -<?+(£)) (1 -$-(£)), 
(E-E ± f 



g± (E) = exp 



2o\ 



(2) 



where all the parameters E±, a±, a±, and 7 are 
iV-dependent; the factors of (I — g±(E)) ensure 
the smoothness of the distribution for E = E±. 
This is essentially a Binder-Landau ansatz cor- 
rected for mixed-phase contributions |§; the re- 
sulting algorithm (MI) works quite well up to 
N = 25, but it becomes very inefficient as N in- 
creases further. 

It turns out that p does not follow the Gaussian 
behavior of Eq. (|J) in the region £L <C E <C 
E + ; if we start from one peak and move towards 
the other, the distribution will eventually follow 
an exponential law. In order to reproduce this 
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Figure 1. The Monte Carlo time evolution of the 
energy, using the M2 algorithm for N = 64. 



behavior, we introduce two new parameters E±, 
E- < E- < E + < E + , and replace g± in Eq. (||) 
with 



g±{E) = exp 



(E-E ± )' 
2a 2 ± 



E < E_ (E > Ej, 



g±(E) = exp(r±E + s±), E > E^ (E < E+), (3) 

where r± and s± are determined by the condition 
that the two branches of g± join smoothly (up 
to the first derivative) for E = E±. The new 
algorithm (M2) works remarkably well up to N = 
64, as shown in Fig. |l]. 

The main drawback of algorithms Ml and M2 
is that the multicanonical parameters must be 
tuned to an ever finer degree with increasing N. 
Thanks to the absence of tunneling in a canoni- 
cal simulation, E± and a± can be estimated accu- 
rately by performing a canonical simulation start- 
ing from a disordered and an ordered configura- 
tion respectively, using the HB algorithm. The 
other parameters can be estimated roughly from 
the simulations at lower N; this estimate needs to 
be refined by performing successive multicanon- 
ical simulations and looking at the resulting en- 
ergy distribution. For N = 64 this required more 
then 10 simulations at moderate statistics (about 
200k sweeps), with a computational effort com- 
parable to the high-statistics simulation with the 
optimized parameters. 
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Table 1 



Summary of high-statistics simulations. 



N 




alg 


stat 


7 




25 


0.3580 


HB 


5M 






25 


0.3574 


Ml 


5M 


5 x 10- 


-3 


3G 


0.3585 


M2 


5M 


1 x 10" 


5 


49 


0.3588 


M2 


4M 


1 x 10" 


10 


64 


0.3595 


M2 


2M 


5 x 10" 


17 



4. Results 

A summary of our high-statistics simulations is 
presented in Table [l]. 

The quality of our ansatz can be judged 
from Fig. |], where the worst case N = 64 is pre- 
sented. The ansatz is not consistent with the data 
within the statistical errors, but it is more then 
accurate enough for the purpose of optimizing the 
multicanonical algorithm. 

It is interesting to notice the exponential fall-off 
of the energy distribution between the two peaks, 
which for N — 64 is followed to great accuracy for 
several orders of magnitude. On the other side of 
the two peaks, the energy distribution falls much 
faster, even faster then the Gaussian behavior as- 
sumed in Eqs. (||) and (||). 

A single multicanonical simulation can be used 
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Figure 2. The (unnormalized) energy distribution 
p, compared with a best fit to Eq. (||). 
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Figure 3. Linear fit in l/N 2 to (it- The solid line 
is a fit excluding the rightmost point (N = 25); 
the dashed line is a fit including all the four points 
plotted. 



to obtain results for a (small) range of j3, using 
the reweighting technique. Our estimator of the 
inverse transition temperature (3 t is the value of [3 
which maximizies the specific heat. Fig. |3| shows 
that 0t is in perfect agreement with the linear 
behavior in 1 /N 2 expected for a first order phase 
transition. Extrapolating to N = oo by a linear 
fit, we obtain 

fit = 0.3596(2). 
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